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1 . Introduction . 

Any reasonable model of a complex communication or other 
service system requires consideration of several interacting 
populations. Consequently, it becomes necessary to name the 
state of each population, and to describe the state of the entire 
system in terms of a vector-valued state variable. With rare 
exceptions very little information can then be derived in a 
direct manner, e.g. by postulating that the modelling process is 
multidimensional birth-death or, more generally, Markov, and 
deriving mathematical expressions for steady-state probabilities, 
etc., as exemplified by Feller [3] I, Chap XVII, and in many 
papers. Exceptions do occur but seem rare, see the cyclic queue 
model of Gordon and Newell [5], and related work by Whittle [10] 
and by Kingman [6] . Consequently, it is tempting to devise 
approximate diffusion models for such processes; previous work wi 
this intent has been done by Schach and McNeill [8], McNeill [7], 
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and the approximation studied by Barbour [2], Our paper, [4], 
indicates how this may be done for a transitory situation. In 
this paper we study a communication system that consists of c 
channels, sought for by arriving messages that balk if they 
encounter an occupied channel. Balking messages enter a distinct 
(retrial) population, or populations, from which they may eventually 
either defect, or again apply for service on one of the channels. 
Previous work in this problem area along classical point-process 
model lines is described by Riordan [9], Our purpose here is to 
indicate how stochastic differential equations may be used to 
write down a model, and how useful information may then be derived. 
In order to illustrate the degree of accuracy achieved we have 
conducted some sample simulations, the results of which generally 
support the approximation method. 
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2. The Model. 



The model structure is depicted in Figure 1. Input to a 

service center arrives according to a Poisson process with 

intensity cX (t) . The service center consists of k distinct 

compartments and a total of c channels or servers. The service 

t h 

process on the i — compartment is Markovian (jjK(t)); i.e. if 
a "customer," here message, is undergoing stage-i service at t 
it terminates within (t,t+dt) with probability y^(t)dt + o (dt) , 

independently of previous process history. Upon completion of the 

t h s t 

i— service function a message proceeds immediately to the (i+1)— . 

The channels are considered separate, and if any compartment of a 
channel is occupied then no other message can enter that channel. 
Consequently, the service center has a capacity of c messages 
at any one time, and a single message has a total service time 
which is the sum of k independent but time-dependent exponentials. 
If y i (t) = y , i=l,...,k then the message service time is 

Gamma (k,y) . In fact, the compartments are introduced in order 
to permit the modelling of non-exponential service. 

We postulate that a message selects a channel uniformly 
at random from among the c possible channels. If that channel 
is busy, the message is denied immediate access to service. Other- 
wise, the message occupies the selected channel until service in 
all k compartments has been completed. The random selection of 
channels and temporary denial of service goes unrecognized in the 
context of classical queueing problems; however, such assumptions 
are quite appropriate for certain communication situations (telephone 
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Figure i 



or air-traffic control, for example) . In such cases the customer 
may not automatically select an empty channel, but must use what- 
ever channel is appropriate, whether busy or not. If channels 
tend to be used equally, the random selection is then appropriate. 

If a customer is not serviced there is no physical queue. 

A telephone caller receiving a busy signal can only hang up, per- 
haps with the intent to call again shortly. An airline pilot or 
submarine attempting to contact a busy traffic controller or shore 
communication facility must try at a later time, possibly to be 
denied service again. We assume that messages or customers who 
are denied service enter a category R. Such customers may retry 
or recall at a later time. These customers reapply for service 
after being in state R for a Markovian (v(t)) - '*' period of time. 
In most cases v(t) would be rather large if the customer must 
receive service (air-traffic control) , but might be small in the 
case of telephone calls. 

In many situations customer discouragement will occur, 

expecially if repeated requests for service are denied. We adopt 

a simple model of customer discouragement by assuming a binomial 

probability 1 - a that a customer upon leaving R decides to 

leave the system and is lost. It is clear that a vast array of 

models can be introduced to represent customer discouragement. 

f 0° 

The most general sort might include states {R.j with R. 

1 i=l 1 

signifying that the customer has been denied service i times. 

The holding time in R^ is Markovian (v^(t)) ^ , and 1 - ou is 
the probability of loss. We discuss only the simple model 
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mentioned earlier with one state R, but the reader may notice 
that the analysis used can easily be extended to handle the most 
general case with states {Ft} 

We wish to describe the behavior of this system by calcu- 
lating system characteristics such as the fraction of customers 
lost, the utilization of the c channels, and the number of 
customers in state R, all as a function of the system parameters 
Standard Markov chain methods can be used on this model owing to 
the Markovian assumptions. Unfortunately, the state of the system 
is a k + 2 dimensional vector. While a Markov chain analysis 
is straightforward, it must be carried out numerically. This 
makes it difficult to assess the influence of system parameters 
on the quantities of interest. Furthermore, it is difficult to 
deduce information about the transient behavior of the system and 
to handle non-stationary transition probabilities using Markov 
chain techniques. For these reasons we adopt the method of diffu- 
sion approximations as an approach. This consists of letting 
c -* +°° and treating the resulting random processes as the sum of 
a deterministic process plus an additive noise (diffusion) process 
This technique allows one to compute the quantities of interest 
as a function of system parameters and to describe the transient 
behavior of the system. If we assume the transitions are station- 
ary, then techniques from the theory of stationary processes, 
especially spectral analysis, can also be applied. While all 
results are exact only in the limit as c -*■ 00 , we present compari 
sons between simulated systems with c = 10 and 20, and diffusion 
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approximations. The comparisons will show that diffusion 
approximations may offer surprising accuracy even for c as 
small as 10. Thus the method described is quite appealing as an 
approximation tool in the present problems, and in many others 
as well. 
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3 . Diffusion Approximation of the System . 



We introduce the following notation: 

Q^(t) = number of messages or customers at service stage i, 
at time t, i = l,...,k, 
k 

Q(t) = £ Q.(t) = number of customers receiving service at 

i=l 1 

time t, 

R(t) = number of customers in state R at time t, 

L(t) = cumulative number of customers lost by time t. 

A customer arriving at the service center at time t selects 
a channel at random, hence with probability Q(t)/c is denied 
immediate service, and with probability l-Q(t)/c begins service. 
It is straightforward to describe this k + 2 dimensional system 
and its transition probabilities over the interval (t,t+dt); 
however, we choose to do this approximately to terms of order dt 
using techniques from the theory of stochastic differential equa- 
tions; see Arnold [1], and also Gaver, Lehoczky and Perlas [4]. 

Using the notation dX(t) = X(t+dt) -X(t) for a stochastic process 
{X (t ) , t £ 0} we express the evolution of our process as follows: 



dQ-^t) = (X (t) c+av (t) R(t) ) (1 - Q(t) /c) dt - ]i 1 ( t) (t) dt 

+ /X ('t)c(l - Q(t)/c) dW x (t) + /ctv (t)R(t) (1 - Q(t)/c) dW R (t) 



- (t) Q x (t) 



dW (t) 



(3.1) 
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dQ i (t) = y i _ 1 (t)Q i _ 1 (t)dt- y i (t)Q i (t)dt+ /y i _ 1 (t)Q._ 1 (t) dW (t) 

y i-l 

“ /y . (t) ( t) dW_ (t) for i = 2,...,k 
1 y i 

dR ( t) = - (l-a+a(l-Q(t)/c) ) v(t)R(t)dt + X(t)c(Q(t)/c)dt 

- / ( 1-a+a ( 1-Q ( t)/c) ) v ( t ) R ( t) dW R (t) + A (t)c (Q(t)/c) dW x (t) 

dL ( t) = (1-a) v (t) R(t) dt + / ( 1-a) v ( t) R ( t ) dW R (t). 

In equations (3.1) the terms dW (t) , dW_(t), and 

U . K 

1 

dW (t) are the "derivatives" of independent standard Wiener 
A 

processes, and as such are usually called Gaussian white noise. 
Other approximations are possible, but are not pursued here. We 
mention, for example, using Poisson white noise where the standard 
Wiener process is replaced by a Poisson process with zero drift. 

A word about the derivation of our equations (3.1) is in 
order. Consider that for the occupancy of the first service 
compartment, (t) : conditional on the values of (t) and 

R(t) the drift or mean change of in time dt is (a) positive 

by the amount (X ( t) c+av ( t) R (t) ) (l-^-^-)dt, where the first 

term represents the expected number of arrivals in (t,t+dt), 
and the other represents the probability of acceptance into ser- 
vice in compartment i = l, while (b) the second, negative, term 
“Uj (t) (t) dt , represents the expected number departing in 

(t,t+dt); hence the difference is the net expected increase in 
. Now for dt small we represent the fluctuating (diffusion) 
component of input by (c) : /X (t) c (1 - Q (t) /c) dW^(t) , where the 
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scale factor is the standard deviation of a Poisson process with 
mean X (t) c (1 - Q (t) /c) and dW^ (t) is W(0,dt), plus (d) a 
corresponding but independent term representing arrivals from R, 
minus (e) another corresponding term representing departure from 
Q^. Although this derivation is heuristic, the rationale is 
simply that in an interval of length dt the various condition- 
ally Poisson components of change are approximately normal, 
owing to the fact that c is presumed large. 

We wish next to view the state of the system 
(Q^ ( t) , . . . , Q k (t) , R ( t) ,L ( t) ) as the sum of a deterministic pro- 
cess plus a noise or diffusion process. To accomplish this we 
introduce the standardized noise variables: 

Q . (t) - cq . (t) k 

= — i = 1 , . . . ,k; X(t) = l X. (t) 

/Z i=l 1 

= R (t) - cr (t) 

✓c 

( 3 . 2 ) 

L (t) - cft(t) 

/z 

k 

= I q ± (t) 
i=l 1 

We then write (Q-^ (t) , . . . ,Q k (t) ,R(t) ,L(t) ) in vector 
fashion as 

c (q^ (t) ,q 2 (t) , . . . ,q k (t) ,r (t) , H (t) ) 

+ c(X 1 (t) ,X 2 (t) ,... ,X k (t) ,Y(t) ,Z(t) ) , 

the first term is the deterministic approximation, the last is 



x i (t) 
Y(t) 
Z (t) 

q (t) 
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the noise process. We substitute definitions (3.2) into (3.1) 
and divide the resulting equations by /c. The results are, to 
terms of order 1 in c, expressible as 



dX 1 (t) = { (1 - q(t) )av(t)Y(t) - (X(t)+av(t)r(t) )X(t) - (t)X x (t) }dt 

+ /X (t) (1 - q(t) ) dW^ + /av (t) r (t) (1 - q ( t) ) dW R - /iJ^TtTq^TtT dW^ 

- (t) - (X ( t) +av (t)r(t)) (1 - q (t) ) + y-^Uq^t) }dt, 

dX i (t) = {y i _ 1 (t)X i _ 1 (t) - y i (t)X i (t) }dt + ✓y i _ 1 (t) q j __ 1 ( t) dW Q 

“ /y • (t)q i (t) dW - /c{q^(t)-y i _ 1 (t)q i _ 1 (t) + y i (t)q i (t) }dt 

w i 

for i = 2, . . . ,k, (3.3) 



dY(t) = {- [ (1-aq (t) ) v (t) ] Y (t) + [ X ( t ) + av ( t ) r ( t) ] X ( t) }dt 

- /(1-aq (t) ) v (t) r (t) dW R + /X ( t) q (t) dW^ 

- /c{r ' (t) + (1-aq (t) ) v (t) r (t) - X ( t) q (t) }dt , 

and 

dZ(t) = (1-a) v (t) Y ( t) dt + / (1-a) v ( t) r (t) dW R 

- SS{V (t) - (1-a) v(t)r (t) }dt. 

We now let c -»• 00 in equations (3.3) . Clearly, in all 
cases the /c term must be identically 0, or else the equations 
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explode. Setting the Sc term to 0 we derive a system of 
ordinary differential equations satisfied by the deterministic 
approximation : 

qj_(t) = (X(t)+av(t)r(t) ) (l-q(t) ) - y 1 (t)q 1 (t) 

q^(t) = P i _ 1 (t)q i _ 1 (t) -P i (t)q i (t) i = 2,...,k (3.4) 

r' (t) = - (1-aq (t) ) v (t) r (t) + X (t) q (t) 

V (t) = (l-o) v(t)r(t) 

It is easy to find q^(t) in terms of q^_^ (t) for 
i = 2,...,k by solving the second equation in (3.4). 

ft ~/ji. (x)dx (s) ds 

q i (t) = j e y i _ 1 (s)q i _ 1 (s) ds + q i (0)e 

but no further explicit results seem attainable without simplifi- 
cation and further parameter specification. Of course numerical 
solution of the differential equations by computer is always 
possible, and may well prove to be the fastest route to useful 
information. 

Steady State Behavior of the System 

Suppose however that we let t 00 and attempt to find a 

steady state solution. As t -* 05 let X (t) , v(t) , y^(t) , 
q^(t), and r(t) converge to X, v, y^ , q^ , and r respec- 
tively. Then q^(t) and r' (t) all converge to 0 and the 
steady state equations become 
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0 = (A+avr) (1-q) - y-j^ 



0 = y . ,q. . - y . q . 

i-l^i-l 1^1 

0 = -(l-aq)vr + Aq 



{,'(<») = (l-a) vr 
k 1 

Letting y = 1/ £ — we find the steady-state values 

i=l y i 



r (co) 



*3 

v ( 1-aq ) ' 



a( cox _ A+y-/(A+y) *-4Aay' 
4 v ' 2ay 



(3.5) 



The negative square root is used for q since 0 £ q £ 1. Fur- 
thermore, the asymptotic loss rate, £'(«>), is given by 
^l ' -aq— . The input rate is A, so the output rate of those 
actually served is A ( > and the asymptotic loss fraction 
(the fraction of customers leaving without receiving service) is 
given by (l-a) q/ (1-aq) . 

We intend to carry out the subsequent analysis under the 
assumption that steady state conditions prevail. In this model 
the deterministic equations are nonlinear, hence there is a ques- 
tion of the stability of the system (in the sense of Liapunov) . 

We must be concerned with the effect of a small perturbation when 
the system is in steady state. If the system is not stable then 
it will diverge from, rather than return to, steady state. The 
stability can be established when A(t), y^(t) and v(t) con- 

verge to A, y^, and v by first linearizing equations (3.4). 
We omit £(t) from the system since it clearly grows without 
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limit and has no steady state value. We assume now that the 
parameters X, }ju , and v are all constants, and express (3.4) 
in matrix form: 





= 


'*1 


+ 


'-(X+y^) -X -X aV’ 




^(t)' 


*2 (t) 




0 




0 
o 

(N 

1 

i—l 






• 




• 




0 y 2 . . . 




• 


• 




• 




• 0 . . 




• 


• 




• 




• 

0 




• 


qk (t) 








O 

1 

• O 

o 




rt 


r’ (t)' 




0 J 




XX X 




r (t) J 



X' (t) = c + AX ( t ) 

**** /vrv 

The system will be stable in the sense of Liapunov if the 
k + 1 eigenvalues of the A matrix have strictly negative real 
parts. The eigenvalues all have strictly negative real parts 
provided vu > 0 and v > 0 for i = l,...,k. We prove this in 
the appendix. 

We now turn to a description of the noise process. We 
assume deterministic equations (3.7) are satisfied. The represent 
tation (3.3) becomes in matrix form 

dU (t ) = A, U ( t ) dt + B.dW. (3.8) 

where U(t) = (X_ L (t) , . . . ,X R (t) ,Y (t) , Z (t) ) • 

W(t) = (w x (t) ,W R (t) ,W Q (t),...,W Q ( t) ) * 

1 k 
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a. 

~t 



■- (y^+X+avr ) 



y 

0 



1 



- (X+avr) 




0 

X+avr 

0 



0 

X+avr 

0 



-(X+avr) av (rq) O-v 

0 0 0 

• • • 

• # 

0 

-y k 0 

X+avr -(l-aq)v 

0 (l-a)v 0 J ' 



and 



B t = rA (1-q) /avr (1-q) -/yq 
0 0 /jJq 



0 

/Xq 

0 



-/ (1-aq) vr 
/(1-aq) vr 



0 

•v/yq 

/yq 

0 



0 

■/yq 

0 

0 



) 



In the above definitions of A^ and we have omitted the 

~t ~t 

argument t in X (t) , y^(t), v(t), q(t), and r(t) for 

notational convenience. As t ■+ <» A^ and B. converge to A 
and B again given by (3.8); however, asymptotically Z(t) is not 
of interest because L(t) -*■ +°°. Furthermore, A. is singular. 

For these reasons we eliminate Z(t) and reduce the dimension to 
k+1. The resulting stochastic differential equation becomes 
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(3.9) 



dV(t) = ^(tjdt + j^d^ 

with V(t) = (X^ (t) , . . . ,X k ( t) ,Y (t) ) ' and given by A with 

the last row and column removed. £ t is given by B^. with the 
last row removed. Again C^_ and D^. converge to C and D as 
t 00 . 

The stochastic process V(t) satisfying (3.9) is a non- 
stationary multivariate Ornstein-Uhlenbeck process. This process 
has been extensively studied, with many results recorded by Arnold 
( [1] , Section 8.2). We shall make use of these results in what 
follows . 

We may integrate (3.9) to find 



V(t) = V(0) + C . V ( t) dt + 
J 0 ^ 



J 0 



D. dW, , 
~t ~t 



(3.10) 



The process V(t) is Gaussian if V(0) is either constant 
or itself Gaussian. This is clear from direct examination of 
(3.10). Suppose V ( t) is either constant or Gaussian. V(t+dt) 
is the sum of V(t) + C J _V(t)dt, which is Gaussian, and another 
Gaussian variable. Thus J7(t+dt) will also be Gaussian and it 
remains to characterize the marginal moments of V(t) as well as 
the covariance structure. 

Let ji t = E(V(t)), J t = E((V(t)-ji t ) (V(t)-Ji t )) ’ be the 
marginal mean and covariance of V(t) . It is shown in Arnold [1] 
that and satisfy first-order differential equations. 

First, the mean vector is described by 



it = £tlit with E t = ^ (0) 



(3.11) 
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is the unique symmetric 



and, second, the covariance matrix £ 
nonnegative definite solution of the matrix differential equation 



it “ £tlt + It£i + BtSi " ith h 



EUViO)-^) mO)-H 0 ) ’) • (3.12) 



Equations (3.9), (3.11), and (3.12) can be formally solved 

with the aid of the fundamental matrix $(t) , that is the matrix 
of solutions of the homogeneous equation $(t) = C^D(t) , $(0) = I,. 

For example, if = C, in steady state, <Ht) = exp (C t) . Using 

4>(t) we find 



t 



V (t) = 4> (t) (V ( 0) + 



J 0 



(* (s) ) _1 D dW ) 

r~ 3 ~S 



Jit = *(t)E(V( 0 )) 



(3.13) 



K(s,t) = E ( (V (s) ~Ji s ) (V ( t) -^ ) ') 



min 



;s,tj 



“ t (s> <Io + 



(*(u) ) 



-1 



D D' ( (4>(u) ) 

~U~U ~ 



-1 



) ’du) $ (t) 



letting s = t, k(s,t) = £ thus providing a formula for . 

We note that in practice it will often be convenient to 

apply a computer routine for solving first-order differential 

equations directly to (3.11) and (3.12). 

Suppose now we assume C. = C and D. = D, that is we 

~t '■“t ~ 

are in steady state. It will be shown in the appendix that all 

k + 1 eigenvalues of C have negative real parts. In fact, the 

matrix C is nearly identical to the matrix A examined earlier. 

f 00 Ct C't 

Then if V(0) ~N(0,7) with Y = e~ DD’e~ dt, V(t) is a 

J 0 

stationary Gaussian process with E(V(t)) = 0, E (V (t) (V(t) ) ’ ) = £ 
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where £ is as defined above or is, equally, the unique nonnegative 
definite solution of the equation 



£i + l£' = -SB’ • 



Furthermore, the covariance function K(s,t) = H(s-t) is given 
by 



H(s-t) 



e C(s t) £ s ^ t ^ 0 

£ e £' (t-s) t ^ s ^ o . 



( 3 . 14 ) 



We summarize our description of the k + 1 dimensional 
queueing system as follows using the diffusion approximation: 



(Q^ ( t) , . . . , (t) ,R (t) ) =r 

C(q 1 (t) , . .. ,q k (t) ,r(t) ) + /C (X 1 (t) , . . . ,X R (t) ,Y (t) ) 

where the first term is given by (3.4) and the second is a multi- 
variate Ornstein-Uhlenbeck process with mean jj^., covariance 
function H(s-t) as described earlier. 



Results for the Single Service Compartment 

We give the exact formulas in the steady state case for 
the special case of k = 1, a single service compartment. In 
this case, the deterministic approximations are still given by 
(3.6). The noise approximation will be a bivariate Ornstein-Uhlen- 
beck process with mean £. The covariance matrix £ will be the 
unique symmetric nonnegative definite solution of the equation 
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-BB' . 



(3.15) 



Al + JV = 



where 



A = 



(- (y+X+avr) (l-q)av 



X + avr - (1-aq) v 



B = 

/V 


r /X (l-q) 


/avr (l-q) 


-/yep 


; BB ' = 


b. b- 
' 1 3> 




/Xq 


-/ (l-aq) Vr 


0 J 




b^ b ' 
^ z 



°22 = Var (Y ( t) ) 

° 12 = Cov (X ( t ) ,Y ( t) ) 



for notational simplicity in what follows. Lastly 



l = ( °ij ) °n = Var(X(t) ) 



Solving (3.15) we find 



°11 = f b i ( |A|+a 22 ) - 2a 12 a 22 b 2 + a 12 b 3 ]/D 

°12 = ^ -a 2l a 22 b l + 2b 2 a ll a 22 ~ b 3 a ll a 12^ //D 
a 22 = ^ a 21 b l ~ 2a ll a 21 b 2 + b 3^ l A l +a 2i-l/ D 
D = 2 ( a n +a 22^ ’ 



(3.16) 



In the last section we present a comparison of the theoretical 
calculations and simulation results in a variety of situations 
with c = 10 and 20. We are approximating the state of 
the system (Q(t),R(t)) by c (q ,r) + /c (X (t) , Y (t) ) . The diffusion 
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approximation allows for more than a description of the marginal 
behavior in that the transient behavior can be characterized and 
the joint behavior at any set of time point, (t^,...,t ) can 
be worked out using the multivariate Ornsfein-Uhlenbeck process. 
However, we do not present numerical results of this kind in 
this paper. 
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4 . The Spectral Matrix . 



We now assume the system is in steady state and adopt the 
steady state diffusion approximation of it. The noise process is 
given by 

dU ( t) = AU(t)dt + BdW(t) . (4.1) 

fv rv /v rv 

The noise process {U(t),t^0} is stationary, hence 
(Q^ (t) , . . . ,R (t) ) will also be stationary. This observation opens 
up the possibility of applying techniques from the theory of 
stationary processes to our queueing process. In particular we 
compute the spectral matrix associated with the noise process. 

We begin with the spectral representations of the stationary 
processes U(t) and W(t) 



U(t) 



iiot . 

e dS (u>) , 



w ( t) = 



ia)t Jr , . . 
e dS w (w) 



(4.2) 



where {S,.(a)),a> e (-<»,<») } and (S n (ut),ii)E(-®,®)} are processes with 

orthogonal increments. The spectral matrices associated with U 

and W, f„ and f r ,, are the matrices of cross spectral densities 
~ ~U ~W 

of the stationary processes and are given by 



£u“ > - E(dS u (o.)dS (0,)), f w ( W ) = EtaS w (a.)dS w (a» ) = 

rv rv #v rw a/ 

We wish to compute f ..(<*>) in terms of A and B. 

rv 

Combining (4.1) and (4.2) by differentiation of (4.2) we 

find 



( iwl-A) dS M 



BdS 7 
~ ~W 



(4.3) 
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Taking transposes in (4.4), multiplying, and taking expectations 
we obtain 

BB ' 

E [ ( iu)I-A) dSydSy ( iwl -A) ] = E [Bdg^dS^B '] = ^ . (4.4) 

rw rv /v /w 

Solving (4.5) for f It = E (dS rl dS TT ) we find 

/V rv /V 

£u (w) = ( ^ ,+ita, i )_1 * (4 - 5) 

We illustrate this computation in the case k = 1, a 



single service 


compartment . 


Here 








r- (y+A+avr) 


(1-q) av' 




(2yq A ■'i 


A = 


A + avr 


- ( 1-aq) v^ 


, BB* = 


• A 2Aq' 



with A = AVq(l-q) - vr/a(l-q) (1-aq) . 

After a few simple matrix inversions and multiplications we 
compute 



II 

^3 


f QQ (a)) 


^3 

a 

< 4-1 




[f RQ (w) 


f RR <“ > 



where 

fQQ (oj) = [2yqto 2 + v (l~aq) (2yqv (l~aq) +Aav (1-q) ) 

+ av (1-q) (Av (1-aq) + 2Aqav (1-q) ] /D 
fR R (w) = [2Aqoo 2 + (A+avr) 2yq (A+avr+A (y+A+avr ) ) 

+ (y+A+avr) ( A (A+avr) + 2Aq (y + A+avr) ] /D (4.6) 
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f QR (u>) = { [Ato 2 + v (1-aq) (2yq (A+avr) + A (y+A+avr) ) 

+ av (1-q) (A (A+avr) + 2Aq (y+A+avr) ) ] 
+ ia) [Av (1-aq) + 2Aqav (1-q) - 2yq (A+avr) 

- A (y+A+avr) ] }/D 

f RQ ( "> = f C.R ( "> 

with 



D = [co 2 + av (1-q) (A+avr) - A (1-aq) (y+A+avr) ] 



2, .2 



+ [v (1-aq) + y + A + avr] a> 



The functions ^qq^' f RR (<v) are the spectral densities 

of Q(t) and R(t) respectively. The real part of f r . r .(o)) 

OR 

is the cospectral density and the imaginary part is the quadrature 
spectral density. The latter two densities provide information 

about the phase behavior of (Q(t),R(t)). Notice that f^tto) 

_2 

and f RR (oj) exhibit tail behavior w and thus fluctuate in a 
high-frequency manner similar to that of the ordinary Ornstein- 
Uhlenbeck process. 
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5. Comparison with Simulation 

The previously quoted results are derived under the supposition 
that c, the number of service channels, becomes indefinitely 
large. Since it will be desirable to apply the approximations in 
case c is finite, we have undertaken to perform several simula- 
tions of particular systems with moderate c-values, and thus to 
provide an empirical check of model accuracy. We shall see that 
our approximations are generally good. 

The simulation simulated is the system (Q (t) ,R(t) ,L (t) ) in 
which k = 1 (the single service compartment case), and X, y , 
and v are constants. By virtue of the Markov nature of the 
system we see that (i) sojourns in states are of independent and 
exponentially distributed duration, and (ii) state changes occur 
in accordance with multinomial Bernoulli trials. Reference to 
( 3.1 ) shows that the state-dependent transition rates are the 
following 



(Q(t) , R ( t ) , L (t ) ) 



(Q(t)+l,R(t) ,L(t) ) 

(Q ( t ) +1 , R ( t ) -1 , L ( t) ) 
(Q(t) -1 ,R (t) ,L (t) +1) 
(Q (t) ,R (t) +1 , L (t) ) 

(Q ( t) ,R(t)-l,L(t) +1) 



Xc [l-Q(t)/c] 

vaR(t) [l-Q(t)/c] 

yQ (t) (5.1) 

XcQ (t) /c 

v(l-o)R(t) 



Hence, given that at time t the system is in state (Q ( t) , R ( t) , L ( t ) ) 
it resides there for independently and exponentially distributed 
times with means equal to the inverse of the sum of the right-hand 
side of (5.1), and then instantaneously jumps to a new, neighboring. 



24 



state with probabilities given by the right side of (5.1) times 
the mean sojourn time in state. Our simulation program is based 
on this scheme. 

The simulation results recorded were the fractions of times 
that the system inhabited each state (Q=i , R=j , L=k) over the 
course of the simulation, the latter chosen to be long. These 
fractions approximate the stationary state probabilities. The 
latter probability estimates were in turn used to compute estimates 
of E[Q («>)], Var[Q(°°)], E[R (<»)], Var[R(°°)], and E[L(»)]. Also, 

the tabulated empirical marginal distributions of simulated Q 
and R were plotted on normal probability paper in order to pro- 

V 

vide a visual check for normality. 

Discussion 

Agreement of the diffusion model with the simulation is, 
apparently, quite acceptable for the cases considered, which by 
design include relatively small numbers of channels. Extensive 
additional simulation results, left unreported here, convey the 
same message. The most noticeable discrepancy occurs in the 
estimates of Var[R(°°)]: that obtained from the analytical 

approximation consistently exceeds the simulation estimate. Attempts 
to show that simulation's failure to reach steady state (starts 
were normally made at Q(0) = R(0) = 0) by starting higher gave 
essentially the same results. The discrepancy remains unexplained. 

In order to aid quick appraisals we have tabulated some key 
quantiles; on the simulation side these are inexact both because 
of the inherent discreteness of the distributions and because of 
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simulation sampling error, and on the diffusion side because of 
the use of the continuous normal approximation to a discrete 
distribution. We have, for example, diffusion approximated 

Q q 95 by E[Q] + 1.65 St. Dev. [Q] , Qq .75 b Y E ^ Q 0.75^ + 

0.68 Std.Dev. [Qq 75 ]/ and Qq by E [Q ] . The selected proba- 

bilities were calculated using a simple continuity correction. 



e.g. 



P{Q * x} 



1 



(X+J-E[Q] ) /Std.Dev. [Q] 

- iz 2 d 2 



/2? J 



— OO 



The latter give at a glance an impression of the Gaussian marginal 
model adequacy, which by and large is good out to the two-sigma 
level. As is to be anticipated, the Gaussian approximation 
degenerates in quality when E[Q] + k Std.Dev. [Q] nears either 
zero or c. 
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Table 1 



Comparison of 


Simulation and 


Diffusion 


Approximation 


II 


II 

c 

II 

<y\ 


a = 0.5, 


c — 10 




Simulation 


Diffusion Approximation 


EQQ (°°) ] 


7.31 




7.34 


Var [Q (°°) ] 


2.00 




2.08 


Std.Dev. [Q (°°) ] 


1.41 




1.44 


E [R (°°) ] 


13.55 




13.54 


Var [R (°°) ] 


16.84 




20.80 


Std.Dev. [R(») ] 


4.10 




4.56 



Q 0.05 


4 


5 


Q 0 . 25 


6 


6 


Q 0.50 


7 


7 


Q 0.75 


8 


8 


o 


9 


10 


0.95 


R 


7 


6 


0.0 5 


R _ _ _ 


10 


10 


0.25 


R . _ . 


13 


14 


0.50 


R _ _ 


16 


17 


0 . 75 


R n „ c 


20 


21 



P{Q (°°) * 9} 


. 956 


.933 


P{Q(°°) £ 7} 


.531 


. 544 


P{Q (°°) £ 5} 


.104 


.102 


P{R(°°) £20} 


.946 


.937 


P{R(°°) i 13} 


.518 


.496 


P{R(°°) ^ 6 } 


.030 


.062 
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Table 2 



Comparison of Simulation and Diffusion Approximation 
X = 5, y = 4, v = 6, a = 0.5, c = 10 





Simulation 


Diffusion Approximation 


E[Q(«) ] 


6.45 


6.49 


Var [Q (°°) ] 


2.36 


2.51 


Std. Dev. [Q ( °°) ] 


1.54 


1.58 


E [R (°°) ] 


8.09 


8.01 


Var [R (°°) ] 


10.20 


12.65 


Std. Dev. [R (°°) ] 


3.19 


3.56 


Q 0.05 


3 


4 


Q 0.25 


5 


5 


Q 0.50 


6 


6 


Q 0. 75 


7 


8 


Q 0.95 


8 


9 


R 0.05 


3 


2 


R 0. 25 


5 


6 


R 0 . 50 


8 


8 


R 0. 75 


10 


10 


R 0 . 95 


13 


14 


P{Q(°°) ^ 9} 


. 986 


.971 


P{Q (°°) £ 6} 


.497 


.504 


P{Q (°°) ^ 3 ) 


.030 


.030 


P{R(°°) £ 14} 


.966 


.965 


P{R(») £ 8} 


.581 


.555 


P (R (°°) £2} 


.021 


. 050 
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Table 3 



Comparison of Simulation and 
X=7, M = 4 , v = 6, 

Simulation 



E [Q (°°) ] 


14 . 62 


Var [Q («>) ] 


4.08 


Std . Dev. [Q (°°) ] 


2.02 


E [R ( 0O ) ] 


27.01 


Var [R (°°) ] 


34.07 


Std. Dev. [R (°°) ] 


5.84 


LD 

O 

O 

a 


11 


Q 0.25 


13 


o 

LD 

O 

a 


14 


Q 0.75 


16 


Q 0. 95 


17 


R _ _ _ 


17 


0.05 


R . _ _ 


22 


0.25 


R. _ . 


26 


0.50 


R _ __ 


30 


0.75 


R 0 . 95 


36 


P (Q (°°) * 19} 


.998 


P {Q (°°) a: 15} 


. 654 


P {Q (°°) * 11} 


.066 


P{R (°°) ^ 40} 


.982 


P { R ( °° ) s; 27} 


. 532 


P {R (°°) ^ 14} 


.026 



Diffusion Approximation 
a = 0.5, c=20 

Diffusion Approximation 

14.69 

4.16 

2.04 

27.08 

41.60 

6.45 

11 

13 

15 

16 
18 
16 
23 
27 
31 
38 

.991 

.655 

.059 

.981 

.526 

.025 
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Table 4 



Comparison of Simulation and Diffusion Approximation 
A = 5, y = 4 , v = 6, a = 0.5, c = 20 



E [Q (°°) ] 



Simulation 



12.89 



Diffusion Approximation 
12.98 



Var [Q (°°) ] 


4.89 


5.02 


Std.Dev. [Q (°°) ] 


2.21 


2.24 


E [R (°°) ] 


15.98 


16.02 


Var [R (°°) ] 


20.58 


25.30 


Std.Dev. [R (°°) ] 


4.54 


5.03 


Q 0 . 05 


9 


9 


Q 0.25 


11 


11 


Q 0 . 50 


12 


13 


Q 0. 75 


14 


15 


Q 0.95 


16 


17 


R 0.05 


8 


8 


R 0 . 25 


12 


13 


R 0. 50 


15 


16 


R 0 . 75 


18 


19 


R 0. 95 


23 


24 


P{Q (°°) 5 : 17} 


.989 


.978 


p{Q(°°) £ 13} 


.596 


.591 


P{Q (°°) & 9} 


.066 


.061 


P {R (°°) £ 24} 


.960 


.954 


p{R(°°) £16} 


.569 


.540 


P{R(°°) S8} 


.037 


.067 
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Appendix 



We prove that the matrices £ from (3.7) and from 
(3.9) have eigenvalues with strictly negative real parts. Both 
matrices have the form 



M 



f-(M 




-c -c 




0 



V 



2 

0 




U 



3 

0 



0 0 0 



-c d' 

0 0 



0 




0 



c c c 



c e 



with c = A, d = av, e = v for A and c = A+avr, d = av(i-q) , 
e = v(l-aq) for In both cases e-d = v(l-a) > 0. 

We wish to solve for the k + 1 roots of the equation 
|M — 01 I = 0. Simple manipulation gives 



|M - 01 1 = 



-(0+y 1 ) 



l 

o 



- (0+U 2 ) 



0 

c 



0 

0 

c 



(0+y k ) 



- (0+e-d) 
0 



0 

- (0+e) 



= (~l) kJ ^ (0+e) it (0+y.) + (-l)^ -1 (0+e-d)D k 

i=l 



k+3 
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where 



D, = c 
k 



— ( 9 +y 2 ) 



0 

0 



(0 + y k _i) 



\-i 



-(0+M k ) 



The determinant D k can be computed recursively as 

k-1 

D k = (9+1, k )D k-i + v 

and these equations can be solved. We find 



D k = 



n ( 9+m • ) - n y. 
i=l i=l 



/ 9 



(A. 1 ) 



Substitution of (A.l) into | M — 0 1, | gives 



0 = (0+e) II (0+y.) + c(9+e-d) 
i=l 1 



n (0+y.) - n y, /0 



i=l 



i=l 1 



(A. 2) 



We wish to show that all k + 1 roots of (2) have strictly 
negative real parts assuming e, e-d, y^,..., y k are all 
strictly positive. 

The k + 1 degree polynomial on the right side of equation 
(A. 2) has strictly positive coefficients, thus by Descartes rule 
only strictly negative real roots are possible. It remains to 
consider the case of complex roots and we assume 0 = a+b^. Since 
complex roots of (A. 2) appear in conjugate pairs we assume b > 0 
without loss in generality. 
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Assuming 0 is complex allows us to rewrite (A. 2) as 



f k y . *\ 

6 ( 9+e ) = c (0+e-d) | IT - lj 



(A. 3) 



We consider the case a ^ 0 or 0 < arg 0 £ in this 



case 



"k y 

" (s4r) 

1=1 *i 



< 1, and it follows that 



Re 



k 


y • i 




t k 


y • 


> 


3tt 


n i 

l i=l 


r 1 1 _ i 


<0 or 


j < arg n ( 
l i=l 


f 1 ) 


_ i 


^0+y. J 1 

l J 


0+y ± 


± 

j 


< X 



Furthermore, 0 < arg (0+e-d) < arg(0+e) < arg(0) £ j, and 
it follows that 



0 < arg(0(0+e)) < arg 



c ( 0+e-d) 



k y . 
1=1 1 



< 2 IT . 



Consequently, 0 cannot be a root of (A. 3) and consequently, (A. 2) 
if a £ 0. We have just proved that any complex root must have 
strictly negative real part, thus all k+1 eigenvalues of M 
have strictly negative real parts. 
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